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ABSTRACT 

We have carried out a hydro dynamical code comparison study of interacting multi- 
phase fluids. The two commonly used techniques of grid and smoothed particle hy- 
drodynamics (SPH) show striking differences in their ability to model processes that 
/\ ' are fundamentally important across many areas of astrophysics. Whilst Eulerian grid 

based methods are able to resolve and treat important dynamical instabilities, such 
as Kelvin-Helmholtz or Ray leigh- Taylor, these processes are poorly or not at all re- 
solved by existing SPH techniques. We show that the reason for this is that SPH, at 
least in its standard implementation, introduces spurious pressure forces on particles 
in regions where there are steep density gradients. This results in a boundary gap of 
the size of the SPH smoothing kernel over which information is not transferred. 

Key words: hydrodynamics - instabilities - turbulence - simulation - astrophysics - 
methods:numerical:SPH - ISMxlouds - galaxies: evolution: formatiomgeneral 



1 INTRODUCTION 

The ability to numerically model interacting fluids is es- 
sential to many areas of astrophysics and other disciplines. 
From the formation of a star and its proto-planetary disk to 
galaxies moving through the intra-cluster medium, dynami- 
cal instabilities such as Kelvin-Helmholz (KH) and Rayleigh- 
Taylor (RT) play a fundamental role in astrophysical struc- 
ture formation. Most popular hydrodynamical methods can 
be divided into two cl asses: techniques following the gas us- 
ing Eulerian grids fe.g. lLanevll998l : lLeveauell998f) and those 
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which follow the Lagrangian motions of gas part icles such 
as ' S moothed Particle Hydrodynamics' (SPH) ([Monaghanl 
1992). Grid based techniques solve the fluid dynamical equa- 
tions by calculating the flux of information through adjacent 
cells, SPH techniques calculate the gas properties on each 
particle by averaging over its nearest neighbours . Due to the 
extensive use of these techniques, it is interesting to carry 
out code comparison studies on well defined problems that 
test their ability to follow the basic gas physics they are 
designed to simulate. 

Our test problem is to follow a dense cold gas cloud 
moving through a low density hot medium. This is specifi- 
cally designed to capture the same physical processes that 
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Figure 1. Illustration of the blob test. The external medium, 
which initially is in pressure equilibrium with the cloud, travels 
with a supersonic velocity creating a bow shock in front of the 
cloud. The post shock flow is subsonic until the smooth flow ac- 
celerates and again obtains supersonic speed on the lateral sides 
of the cloud. 



occur during the formation and evolution of astrophysi- 
cal structures. We will also study the shearing motion of 
two fluids of different densities to elucidate the problems 
that we find with this test. Similar configurations, includ- 
ing Sh ockwave interaction with clouds, have been studied 
by e.g . lMurrav et alJ (ll99^.lKlein et al.l (ll994MVietri et all 
(|l997ftnMori Burkertl (|20Q(]| ). However we are not aware 
of a direct comparison between simulation methods in this 
context. Differences were found in the literature between 
different studies of the same problem. For example, SPH 
studies of galaxy- cluster interactions by Abad fet alJ (1999) 
found that only half the inter-stellar medium was removed 
from the galaxy. Using a grid based calculation of the same 
initial conditions, Quil is et alJ (|2000) found that all the gas 
could be removed and attributed this to the high resolution 
shock capturing ability of their Eulerian code. 



2 THE BLOB TEST 

A schematic view of the blob test problem can be seen in 
Fig. Q A spherical cloud of gas is placed in a wind tunnel 
with periodic boundary conditions. The ambient medium is 
ten times hotter and ten times less dense than the cloud so 
that it is in pressure equilibrium with the latter. We will 
refer to this initial density contrast between the cloud and 
the medium as Xini- All of the gas is atomic hydrogen with 
molecular weight m mo i = 1.0 and an adiabatic index 7 = 
5/3. 

This setup will investigate how different simulation 
codes handle typical astrophysical processes important for 
multi-density and multi-phase systems, such as ram-pressure 
stripping and fragmentation through KH and RT instabili- 
ties. 



Figure 2. We plot the Mach number of the flow directly down- 
stream of the shock on the symmetry axis of the cloud. The flow 
speed increases due to the weakened shock strength up to £ S onic 
where the relative motion of the cloud and wind turns subsonic. 



3 ANALYTICAL EXPECTATIONS 

Although the nonlinear stages of the KH and RT instabili- 
ties cannot be fully described analytically, we can still use 
analytic arguments to estimate the characteristic disruption 
timescale for the cloud. 

In order to specify our problem we characterize the ex- 
ternal medium with a sound speed c s and assign it an initial 
velocity v — Mc s with Mach number M — 2.7. Further- 
more, we place the cloud initially at rest in the computa- 
tional domain. Since the wind is supersonic, a bow shock 
will form in front of the cloud with the post shock proper- 
ties given by the Rankine-Hugonoit shock jump conditions. 
Because the cloud is accelerated by the wind, we will from 
now on perform all of our calculations in the rest frame of 
the bow shock, referring to pre-shock quantities with the sub- 
script 1 and post- shock with 2. The shock co nditions fo r the 
density, velocity and Mach number are (e.g.|ShJl992) 



p2 

Pi 



Vl 
V2 



(7 + l) + (7-l)(X?-l) 



M 2 2 = 



2 + (7" 



l)Mi 



2jMl - (7 - 1) 



(1) 



(2) 



Formally we would take the obliqueness of the bow shock 
into account but for simplicity we will only consider the 
flow that enters at the symmetry axis of the cloud. 

The cloud acceleration can be approximated by consid- 
ering the maximum area that can gain momentum from the 
ambient flow. This implies that all gas in a cylinder in front 
of the cloud transfers momentum leading to an acceleration 



del ~ vi 



(3) 



Integrating this equation leads us to the evolution of the 
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pre-shock velocity 



«i(t) 



(t + Z/Vext)' 



(4) 



where / is a characteristic length given by / = M c \/27rR 2 l p e xt- 
By using Eq. 0]to calculate the pre-shock Mach number to- 
gether with Eq. 0we can obtain a qualitative understanding 
of the post-shock velocity. This velocity is crucial for the sta- 
bility of the cloud surface and, as we will show in section 13.11 
for the destruction of the cloud itself. The evolution of the 
post-shock Mach number M2 is given by 
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for t <C tsonic 
for t > tsonic 



(5) 



Here tsonic is the time at which Mi — M2 — 1 and the shock 
disapperars. After this point, gas freely streams towards the 
cloud and the Mach number decreases only due to the con- 
tinued acceleration. Notice that for t < t SO nic, M2 < 1, even 
for Mi = vi/cs — > 00. This means that behind the shock, 
the flow will always be subsonic and we expect instabilities 
to grow there. For t — > 00, M2 — > and the cloud will even- 
tually be co-moving with the background flow. The evolution 
of the postshock Mach number is shown in Fig.|2]in terms 
of the so-called "crushing time" defined as, in our notation, 




Figure 3. The time dependence of the growth rates of KH (solid, 
blue lines) and RT (dashed, red lines) instabilities. The lines rep- 
resent different sizes of perturbation wavelengths: R c \ (thick), 
R c i/2 (middle) and R c i/3 (thin). 



1/2 



Vl 
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where x is the density contrast between the cloud and the 
external medium. This is a natural timescale supersonic 
cloud evolution. We will naively use % = Xmi = 10 and 
vi = v ex t, representing our initial condition. During the in- 
terval of T cr a bowshock is formed and the shocked gas will 
form a smooth flow around the cloud, reaching supersonic 
speed at the points indicated in Fig. [I] Beyond this region 
we expect to see a turbulent boundary layer forming which 
transports material off the surface. The cloud will compress 
along the line of motion due to an internal shock wave gener- 
ated by the external gas. From Bernoulli's theorem we know 
that the pressure is low on the lateral sides which causes an 
overspilling of the cloud due to the high inner pressure of 
the compressed cloud. This causes mass loss irrespective of 
any instability. 



3.1 The Kelvin-Helmholtz instability 

Kelvin-Helmoltz instabilities (KHI) occur when velocity 
shear is present at the interface between two fluids. The 
importance of the KHI, in the context of ga s cloud stabi l- 
ity, has been s tudied by many authors e .g. |Nulsen] (^982), 
IMurrav et all (I1993T) . IVietri et all (Il99/t) . iMori fc Burked 
fcOOfl . 

Neglecting gravity, the dispersion relation of the KHI, 
in the notation of our setup, for an incompressible fluid is 
(IChandrasekharlll96lh 
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where k is the wavenumber of the instability and the last 
approximation holds for x ^> 1. The characteristic growth 



time for the KHI is then 
2tt 
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By naively using the post-shock quantities of Eq.QJand our 
choice of cloud parameters, we can calculate an approximate 
time dependence of the KH instability, which is shown in 
Fig. |3| (blue, solid lines) for perturbations of size R c \ (thick), 
R c \/2 (middle) and R c \/S (thin). Small scale instabilities 
grow faster due to the tkh ~ k -1 relation. The first modes 
to grow are the shortest. Their growth will act to widen the 
interface between the shearing layers, hence dampening the 
growth of modes smal ler than the thickness of the interface 
( Cha ndrasekharlll961r) . The fastest growing modes are now 
those that are equal to the thickness of the interface. As 
this process continues, the mode responsible for the cloud 
destruction is that which is comparable to the size of th e 
cloud itself: k c \ ~ 2tt / 'R c i (lNulsenlll98llMurrav et al.lll993T) . 

The instability growth time is always larger than the 
cloud crushing time. The horizontal line at r — 1.6 r cr in 
Fig. |3 indicates roughly the time at which the k c \ KH mode 
should have grown fully. We will from now on refer to this 
time as tkh. 

Note that cloud compressibility can be taken into 
account when calculating the KH growth time (see 
but was omitted for simplicity. Also 



^ikhl inin ■ 

note that in certain more physically motivated situations 
with external gravitational fields, self gravity, physical vis- 
cosity, magnetic fields, radiation etc., the KHI is modi- 
fied and is damped in many cas es (e. g. IMurrav et al.lll993l : 
IVietri et alJll997t iMiniati et alilTooS iGregori et al.ll200(tl . 

3.2 The Rayleigh-Taylor instability 

Ray leigh- Taylor instabilities (RTI) occur when a denser 
fluid is accelerated by a less dense fluid, or when a heav- 
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ier fluid displaces a lighter fluid. The cloud is accelerated 
with respect to the background and we expect RT insta- 
biliti es to devek>p^Th e dispersion relation for the RTI is 

\w 2 \=k'a( pcl ~ pext ) «Jfe'a, (9) 

where the last approximation is valid for % ^> 1. The KHI, 
which results from shearing flows, has a 2D geometry, and 
can be descibed by as single wavevector k. By constrast, 
the RTI necessarily has a 3D geometry and must be de- 
scribed by a vector wavelength, k' — (ki, fe), of magnitude: 
k' — \Jk\ +k\. The acceleration on the surface can be as- 
sumed to be a = ea c \, where a c \ is given by Eq.|3]and e is an 
efficiency factor. Note that it is very difficult to analytically 
determine the efficiency of the momentum transfer from the 
external medium onto the cloud. By using e = 1 we will get 
a lower limit on trt- 

Fig. |3| shows, for our choice of parameters, the charac- 
teristic growth times for RT instabilities (red, dashed lines) 
of size R c \ (thick), R c \/2 (middle) and R c \/3 (thin), demon- 
strating that tkh < trt for large instabilities. The largest 
mode grows very slowly and is probably not important in 
this type of problem. However, we expect that a fast grow- 
ing small-scale RT instability should develop on the cloud 
front, especially on the axis of symmetry as the flow rams 
into the stagnation point. Complicated mixtures of KHI and 
RTI during later evolution is also expected until the cloud 
becomes fully comoving with the flow. 



4 NUMERICAL SIMULATIONS 

Our numerical simulations solve the Euler equations which 
neglect physical viscosity and radiative processes; we assume 
a perfect gas equation of state. Away from shocks, the evo- 
lution is strictly adiabatic. This means the gas can only 
undergo heating and cooling by adiabatic compression or 
expansion, or by irreversible heating at shocks. In order to 
isolate the differences in hydrodynamic solvers, we neglect 
the self gravity of the gas. 



bilities of the cloud. Formally this can be seen as a triggering 
of small scale RT and Richtmyer-Meshkov instabilities 1 . 

This particle setup is used as IC for the SPH simu- 
lations. The ICs for the grid simulations are obtained by 
smoothing the gas quantities (density, temperature and ve- 
locities), using the standard SPH kernel with 32 nearest 
neighbours, and then mapping these on to a uniform grid. 
In this way the noise introduced by using discrete particles 
in the SPH simulations is also present in the grid IC. As 
we will argue below, the key parameters to study are those 
connected with the resolution and strength of artificial vis- 
cosity therefore our parameter space studies will focus on 
the effect of these. 



4.2 The codes 

The simulation was carried out with about a dozen differ- 
ent independent simulation codes. Since all the grid codes 
gave consistent results, and similar for the SPH codes, we 
shall just present the detailed analysis of a selection of these 
codes which are summarized in Table Q Here we give a brief 
description of these codes and the methods used for solving 
the hydrodynamical equations: 



4.2.1 ART (AMR) 

ART (Adapti ve Refinement Tree) is a iV-body+ gasdynamics 
AMR code dKravtsovl ITooS iKravtsov et all I2002T) . The 
ART c ode uses second-order shock-capturing Godunov-type 
solver (IColella fc Glad 119851 to compute numerical fluxes 
of gas variables through each cell interface, with "left" and 
"righ t" states estima ted using piecewise linear reconstruc- 
tion (jvan Leerlll979l) . This is a monotone method that is 
known to provide good results for a variety of flow regimes 
and resolves shocks within ^1 — 2 cells. A small amount of 
dissipation in t he form of artificial d iffusion is added to nu- 
merical fluxes (|Colella &; Woodw ard 1984), as is customary 
in the shock-capturing codes. The details of the flux eval- 
uation and summation on mesh interfaces can be found in 
iKhokhlovl (|l998F) . In the simulations presented in this paper, 
a new distributed MPI version of the ART code developed 
by Douglas Rudd and Andrey Kravtsov was used. 



4.1 Initial conditions 

The initial conditions (IC) for the blob test are set up in 
the following way: we use a periodic simulation box of size, 
in units of the cloud radius R c \, {L x , L y , L z } = {10, 10,40} 
where we center the cloud at {x, y, z} = {5, 5, 5}. The ICs 
are generated by randomly placing equal mass particles to 
obtain the the correct densities and cloud radius. Using an 
SPH code, the system is evolved and allowed to relax to ob- 
tain pressure equilibrium. Velocities are added to particles 
that are in the surrounding low dense ambient medium. One 
could smoothly increase the velocities to be more faithful to 
astrophysical situations, but this more violent start together 
with particle noise serves as the initial seed for surface insta- 



4.2.2 CHARM (AMR) 

CHARM is an iV-body+gasdynamics, AMR code, based on 
the CHOMBO-AMR library, employing a higher order Go- 
dunov 's method for the solutio n of the hydrodynamic equa- 
tions (Mini ati fc Colellal |2006 ) . Here a piecewice linear re- 
construction scheme with Van Leer's limiter and a nonlinear 
Riemann solver were used, resulting in a second order accu- 
rate method in both space and time. CHARM was used to 
test the influence of IC on the cloud evolution. 



1 The Richtmyer-Meshkov instability occurs when a contact dis- 
continuity gets shocked or rapidly accelerated. This generates vor- 
ticity and structures similar to those of RT (e. g. llnogamovl ir999 ) . 



Simulating fluids using SPH and grid techniques 5 



4.2.3 Enzo (AMR) 

Enzo is an Eulerian AMR hybrid code (hydro + TV- 
body) code that was originally written by Greg Bryan 
and Michael Norman at the National Center for Su- 
per computing Applicat ions at the University of Illinois 
([Brvan fc jj^rnian| |l997|l. Enzo uses the Piecewise parabolic 
method (|Colella Sz Woodwardl 17984) for solving fluid equa- 
tions. PPM is a higher order accurate version of Godunov's 
method with an accurate piecewise parabolic interpolation 
and a non-linear Riemann solver for shock conditions. The 
method is third order accurate in space and second order 
in time. This together with the Riemann solver results in a 
very accurate shock treatment compared to the SPH codes 
where artificial viscosity is used. 



4.2.4 FLASH (AMR) 

FLASH is an AMR hybrid code (hydro + TV-body) de- 
veloped by the AS C Center at the University of Chicago 
JFrvxell et al.l2 000V The PPM hydrodynamical solver is for- 
mally accurate to second order in both space and time but 
performs the most critical steps to third- or fourth-order 
accuracy. For the simulations performed in this paper we 
have used FLASH version 2.3 using AMR with maximum 
refinement up to the resolutions indicated in Tabled 



4.2.5 Gasoline (SPH) 

Gasoline is a par allel Tree + SPH code, described in 
Wad slev et alJ (|2004l) . The code is an exten sion to the N - 
Body gravity code PKDGRAV developed bv lStadell (l200lh . 
Gasoline uses artificial viscosity (AV) to resolve sh ocks and 
has a n implementation of the shear reduc ed version (jBalsaral 
1995) of the standard (jMonaghanlll992T ) artificial viscosity. 
Gasoline solves the energy equation using the asymmetric 
form and conserves entropy closely. It uses the standard 
spline form smoothing kernel with compact support for the 
softening of the gravitational and SPH quantities. 

The AV is implemented by solving a momentum equa- 
tion of the form 



dvi 
~dt 



(10) 



where Pj is pressure, v« velocity and the AV term liij is 
given by, 



where r 



4(c i +c J )/^ J +/3/xf ■ 
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for Vij • r i: 
otherwise, 

h(Vij • Ti 
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Ti — Tj, Vij = Vi — Vj and Cj is the sound 
speed, a and are the coefficients used for setting the vis- 
cosity strength, and are essential for capturing shocks and 
preventing particle interpenetration. Note that the viscosity 
term vanishes for non-approaching particles. The commonly 
used values in the literat ure is a = 1 and = 2 which 
originally was proposed bv lLattanzio et ail (|l986l) using Sod 
shock tube tests. Later we will carry out experiments with 
different values of a and 0. 



Table 1. Simulation details. Enzo and ART use the static grids 
indicated in the table while the CHARM and FLASH simulations 
have been run using AMR up to the indicated resolution. All grid 
simulations are started from the 256,256,1024 initial conditions 
except for the analytically started CHARM run. 



nParticles/Grid size 


AV 


Name 


ART, static 
64,64,256 
128,128,512 
9^fi 9^fi 1 094 


no AV 
no AV 
no AV 


ART_64 
ART.128 
ART 9^(S 


CHARM, AMR 

512,512,2048 


no AV 


CHARM.512 


Enzo, static 

64,64,256 

128,128,512 

256,256,1024 

256,256,1024 


no AV 
no AV 
no AV 
ZEUS, a = 1, = 2 


Enzo_64 
Enzo_128 
Enzo_256 
Enzo_ZEUS 


FLASH, AMR 

64,64,256 

128,128,512 

256,256,1024 


no AV 
no AV 
no AV 


FLASH.64 

FLASH.128 

FLASH.256 


Gadget-2 

10 7 


a = 0.8 


Gad_10m 


Gasoline 

10 6 
10 7 
10 7 
10 7 
10 7 
10 7 


a = 1, = 2 
a = 1, = 2 
a = 0, = 2.0 
a = 0, = 0.5 
a = 0, = 0.1 
Balsara, a = 1, = 2 


Gas_lm 

Gas_10m 

Gas_10mAVl 

Gas_10mAV2 

Gas_10mAV3 

Gas_Bals 



4.2.6 GADGET-2 (SPH) 

The TreeSPH cod e G ADG ET-2 

(jSpringel. Yoshida Whitel l200lL ISoringeJ 1200,^ is 
the updated version of the GADGET- 1. The code is similar 
in character to Gasoline but uses an entropy conserving 
formulation of SPH. This means that the thermodynamic 
state of each fluid element in GADGET-2 is defined through 
the specific entropy and not the specific thermal energy. 
GADGET-2 uses a somewhat different formulation of 
artificial viscosity than Gasoline. The viscosity term in 
EqQSlis here formulated as 



Pij 



(13) 



where v s -] 9 



— Ci ~\~ Cj 



3wij is the so called signal velocity. 



Here Wij = • Yij/\vij\ is the relative velocity projected 
onto the separation vector provided particles approach each 
other, otherwise the term vanishes just like in Gasoline. 



5 RESULTS OF THE SIMULATIONS 

Fig. 4 shows central density slices of Gasoline (Gas_10m), 
GADGET-2 (Gad_10m), Enzo (Enzo_256), FLASH 
(FLASH.256) and ART-Hydro (ART.256). These are the 
high resolution simulations with the default standard 
settings. 
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Figure 5. A thin central slice of the SPH particles of Gas_10m at t = 0.75rKH(left panel) and t = 1.5 tkh (right panel). The density 
ranges from high (blue) to low (white) and the magnitude of the velocity vectors are normalized to the reference frame of the centre of 
the cloud. We clearly see the effect of the cloud stretching due to the lateral Bernoulli zones and the formation of downstream vorticity. 



The simulations of the two SPH codes, Gas_10m and 
Gad_10m, show a very similar evolution. As expected, a 
detached bow shock forms directly in front of the cloud. 
An internal shockwave forms within the cloud compress- 
ing it. The post shock flow encompasses the cloud, creating 
Bernouilli zones on the top and bottom with lower pres- 
sure. This causes the cloud to become elongated as well as 
compressed along the z-axis and we see gas being ablated, 
i.e. stripped through the induced pressure differences, from 
the top and bottom edges. Gas stripping slowly progresses 
and the cl oud's shape does not change significantly for a 
long time ([Doroshkevich &: Zeldovichl ll98lh . Fig. shows 
the particles in a thin slice centred on the cloud. The veloc- 
ity vectors of each particle are plotted in a reference frame 
centred on the cloud. The colours indicate the gas density. 
Behind the edges of the cloud we see a vortex created due to 
the shearing motion of the ambient medium which creates 
a low pressure region behind the cloud. Initially, the cloud 
evolution is similar in the grid simulations. It is compressed 
and elongated and gas is removed from the trailing edges 
where the vortex has created a vacuum behind the cloud. 
Some of the ambient medium turns around and falls onto 
the backside of the cloud. However, the late cloud evolution 
is quite different in these simulations. Early on we observe 
surface perturbations on the front of the cloud, probably 
originating from the way the ICs are setup (see argument 
in sect . 14.1 J) . A complicated mixture of KHIs and RTIs are 
developing on the cloud front which, due to subsequent com- 
pression and lateral expansion, becomes even more KH and 
RT unstable. By t ~ tkh, large scale KHIs have developed 



and the cloud starts to fragment. Further instabilities and 
turbulence mixes the smaller clumps of gas into the ambi- 
ent medium. All grid simulations show basically the same 
cloud destruction time. We also note that Eulerian (shock 
capturing) methods effectively localise shocks to a few grid 
cells compared to the smoothed out shocks in the SPH sim- 
ulations resulting from AV shock capturing schemes. 



In Fig. El we show the remaining cloud mass fractions 
as a function of time for the Enzo and Gasoline simulations. 
These are representative of grid and SPH methods. We de- 
fine the cloud as being any gas that satisfies T < 0.9T ex t 
and p > 0.64 p c i. It is of course possible to construct more 
elaborate criteria but these select the gas that visually is a 
part of the cloud. The figure shows that both techniques give 
a similar mass loss up to ~ tkh • Before this time the gas loss 
is mainly due to ablation into the low pressure zone created 
behind the cloud. As soon as we pass tkh for large scale 
KHI the SPH and grid methods diverge. In the grid simula- 
tion, the cloud quickly disrupts and diffuses into the ambient 
medium, while the SPH simulation only shows continuing 
stripping. After t = 2.5tkh, no gas in the grid simulation 
can satisfy our criteria while the SPH simulation still shows 
a mass fraction of ~ 40%. This shows us that the vortex 
shedding through the Bernoulli zones is the most important 
mechanism for mass-loss at t < tkh in both methods. After 
this time dynamical instabilities dominate the grid mass- 
loss. 
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Figure 7. Resolution study for Enzo and Gasoline. The panels show density slices of, from top to bottom, Enzo_64, Enzo_128, Enzo_256, 
Gas_lm and Gas_10m for t = 0.25,0.75,1.5 and 2.25 tkh- We see that resolution changes the phase of the instabilities in the grid 
simulations while the destruction time is the same. Higher resolution also shows less diffusion and better resolves small scale fragments. 
The Gasoline runs are not able to resolve small scale instabilities at all. 





t/r m 

Figure 6. The evolution of the cloud mass fraction. In the SPH 
simulation (solid, red), the cloud slowly loses mass to the ambient 
medium and has not been completely mixed even after 5 tkh- The 
grid simulation (dashed, blue) follows the SPH up to the time at 
which the KH instability causes it to rapidly fragment and mix. 

5.1 Resolution dependence 

It is difficult to do a direct translation between grid and 
SPH resolution. The maximum allowed resolution is a fixed 
grid of size 256x256x1024 in the grid runs and 10 7 particles 
in the SPH runs. This means that there is almost a factor 
7 more cells compared to particles. On the other hand, cells 
are uniformly distributed in space and only « 21378 cells 
cover the cloud, assuming an IC setup. This should be com- 
pared to the 100, 000 particles constituting the cloud in the 
high resolution SPH run. A comparison like this is still not 
straightforward due to the fact that SPH uses particles as 
non-independant resolution elements. This means that each 
particle is not a carrier of information without neighbours 
to smooth over, and the effective number of resolution ele- 
ments is more or less set by the kernel shape and number of 
neighbours to smooth over. 

Resolution affects the convergence of hydrodynamical 
simulations. A cut-off is always introduced on the scale of 
the spatial resolution below which instabilities can not be 
resolved. This often serves as a source of numerical viscosity. 

For most of the codes used in the comparison, we have 
varied the resolution in order to obtain an understanding of 
how this changes the cloud morphology, mass loss and frag- 
mentation time (see TableCJ. Fig. shows the outcome of, 
from top to bottom, Enzo_64, Enzo_128, Enzo_256, Gas_lm, 
and Gas_10m. In the grid simulations we conclude that, 
while the compression and elongation of the cloud is rel- 
atively similar, the detailed way the cloud frag ments is reso- 
lution dependent, owing to the differences in IC (I Jones et alJ 
1996). In Enzo_64, a phase of the KHI centered on the clouds 
symmetry axis is dominant. This phase is less pronounces 
as resolution is increased (Enzo_128) and it is nowhere to be 
seen in In Enzo_256. Going to higher resolution we see more 
and more small scale instabilities developing which enhance 
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Figure 8. The evolution of the cloud mass fraction for different 
resolutions. As the resolution of the grid simulations is increased 
from 64 to 128 to 256 cells across the wind tunnel, the amount of 
mass increases a little but converges. Increasing the resolution of 
the SPH simulations does not decrease the amount of mass lost, 
rather the opposite, perhaps due to the momentum transfer due 
to massive particles acting like '"bullets'". 

mixing of the cloud material with the background flow. Nu- 
merical diffusion is stronger in low resolution simulations 
which is why parts of the cloud survive longer in the higher 
resolution runs. 

The different SPH simulations are qualitatively very 
similar. Instabilities can not be resolved in Gas_lm nor in 
Gas_10m. However, we note a weak large scale RT instabil- 
ity on the cloud front at t = 2.25 tkh in Gas_10m, which is 
absent in Gas_lm. 

The general description above is again quantified by 
studying the cloud mass fraction at each timestep, see Fig-IHI 
In this plot we have also added an extra low resolution SPH 
simulation using only 100 000 particles. The grid simulations 
show a clear trend of dissolving the cloud quickly after ~ 
tkh regardless of resolution while the SPH simulations only 
show a steady mass loss due to the material ablated into 
the trailing vacuum. Decreasing the SPH resolution causes 
the mass fraction to rise above the initial value during the 
initial phase and mass is lost more rapidly for t > tkh- The 
latter effect is most probably due to the increased mass of 
each particle, causing each particle interaction to transfer 
momentum in a more violent, "bullet-like" fashion. 

5.2 Initial Seeds 

As partly shown in the previous test, the development of 
the instabilities, particularly during the nonlinear stages, is 
sensitive to the exact definition of the initial conditions. This 
is because they set the seed perturbations out of which the 
instabilities grow. However, while the mixing of the cloud 
material with the background medium is affected by small 
scale motions that arise from the small unstable scales, the 
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Figure 9. Evolution of the cloud with 'analytic' initial conditions. Each frame shows a density slice through the cloud center at times 
t = 0.24,0.9, 1.7 and 2.5 tkh with densities varying from low (red) to high (blue). 



cloud disruption is mostly the result of the development of 
the large scale perturbations. 

As an example of this in Fig. |5|we show the evolution 
of the cloud-wind interaction but with initial conditions set 
directly from the analytic definition. Thus in this case the 
initial conditions are free of noise and are purely symmetric. 
A base grid of (32 x 128 x 128) was used with two additional 
levels of refinement with refinement ratio of 4 placed dynam- 
ically in regions where the relative change in density, Ap/p 
exceeded 20%. This corresponds to an effective resolution of 
512 x 512 x 2048 in the finest grids, which reduces the level 
of perturbation with respect to the previous cases. 

As shown in panel B of Fig. |5|the most destructive mode 
has a different phase than in the cases illustrated above for 
the corresponding grid based codes. However, as in the pre- 
viuos cases, by t = 2.5 tkh (panel D) the cloud has been 
completely reduced to debris by the instabilities. This shows 
that despite differences in the appearence of the cloud gas 
distribution its fundamental fate of disruption and subse- 
quent mixing on a timescale of a few tkh is independent of 
the specific definition of the initial conditions. 



6 WHY SO DIFFERENT? 

What is the reason for the observed discrepancies be- 
tween simulations carried out using SPH and grid-based 
techniques? Differences between SPH and grid-based 
results have been discussed before in the literature 
JPearce et allll99flt iThacker et al1l2000l: iRitchie fc Thomas! 
20011: iTittlev et al.1 l200ll : ISpringel fc Hernauistl 120021: 
Marri &; White! I2003F) in different contexts to this study. 
While artificial viscosity is the most obvious focus for 
criticism of SPH it is not the main reason for the differences 
observed in this test. We will show this in section IETT1 before 
focusing on the almost complete suppression of KH (and 
RT) instabilities in SPH simulations of this test and present 
an explanation of why this occurs. 



6.1 Artificial viscosity 

The artificial viscosity (3 parameter in Ea. llll is necessary for 
shock capturing and is required for SPH to work properly 
in supersonic regimes. The a parameter has a less obvious 
meaning and the classical a = 1.0 se tting is most proba - 
bly unphysical. It can be argued (e.g. IWatkins et al.l [1996") 
that a can roughly be interpreted as a Navier- Stokes shear 



plus bulk viscosity, even though the AV is only sensetive to 
flow properties such as interparticle travelling. Bulk viscos- 
ity is normally not important in fluid dynam ics, except i n 
the theory of attenuation of sound waves fe.g. lFabeil ll995 > ). 
In numerical simulation its inclusion is for the most part to 
dampen the so called post-shock ringing. 

The inclusion of artificial viscosity leads us to one of 
the first possibilities for the observed discrepancy: We are 
not solving the same hydro- dynamical equations in the dif- 
ferent codes. By adding AV we are solving some kind of 
Navier- Stokes equation when we actually want to compare 
the solutions to the grid codes that, in this sense, are closer 
to the Euler equations. Note however that there is always 
numerical viscosity due to resolution and truncation in all 
simulation methods. 

Viscosity has two major effects on the processes we want 
to capture in this test: 

(i) Dampening of small scale velocity perturbations and 
random velocities 

(ii) Diffusion of post shock vorticity and thus smearing of 
turbulence 

The effect of (i) will enter as a stabilizing factor for the 
growth of instabilities. Physical kinematic viscosity, is, sets 
a cut off f or the size of the smallest eddies in turbulence 
below which turbulent motion is diffused. 
The effect of (ii) follows from the first one and is ob vious 
from inspection of the vorticity transport equation (e.g. |Shul 

ET 



— + Vx (w x v) 



VP x V ( ij +zyV 2 cj, (14) 



where uj = V x v is the vorticity. The two terms on the right 
hand side can create or diffuse vorticity. The first term is the 
baroclinic term which is non vanishing if we have non-aligned 
pressure and density gradients. This is the case in oblique 
shocks like in the bow shock of our cloud simulation. The 
second term is responsible for diffusing vorticity in space i.e. 
taking local vorticity and spreading it into the general flow. 
This means that as soon as we have viscosity, we will dampen 
vorticity. Especially important is the vorticity generated in 
the post shock flow, which should act to destabilize the cloud 
together with the surface instabilities. 

A study on how ar tificial viscos i ty dam pens small scale 
vorticity was made bv iDolag et alJ (2005). By using a low 
viscosity formulation of SPH they find higher levels of turbu- 
lent gas motions in the ICM and noted that shocked clouds 
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Figure 10. Viscosity study for Gasoline. The panels show density slices of, from top to bottom, Gas_10m, Gas_Bals, Gas_10mAVl, 
Gas_10mAV2 and Gas_10mAV3 for t = 0.25,0.75,1.5 and 2.25 tkh- We can see how reducing shear viscosity and removing the bulk 
viscosity renders very similar results; the cloud destabilizes to a higher degree. By reducing the shock capturing viscosity the cloud 
destabilizes even further, most probably to an unphysical solution in the lower setting. The artificial post shock ringings also gets more 
pronounced, as expected for lower viscosity settings. 



12 Oscar Agertz et al. 



Table 2. Performed KH runs 



Resolution 


X 


<W^shear 




IC 


Name 


Enzo 












{256,256,8} 


8.0 


1/80 


1.70 


lattice 


GRID1 


{256,256,8} 


10.0 


1/40 


1.86 


poisson 


GRID2 


{256, 256, 8} 


10.0 


1/40 


1.86 


glass 


GRID3 


Gasoline 












900 k part 


8.0 


1/80 


1.70 


lattice 


SPH1 


1.1 M part 


10.0 


1/40 


1.86 


poisson 


SPH2 


1.1 M part 


10.0 


1/40 


1.86 


glass 


SPH3 



tend to be unstable at earlier times. However, by looking 
at their Figure 3 we note that the overall difference in the 
cloud evolution is small. 

In order to understand the effect of artificial viscos- 
ity in our cloud-wind test we have performed three sim- 
ulations with modified setting of the viscosity coefficients. 
These are Gas_10mAVl, Gas_10mAV2 and Gas_10AV3, see 
Table Qfor viscosity values. A simulation using the Balsara 
switch but with the standard (a = 1.0, f3 = 2.0) was also 
performed. Fig. 1101 shows the outcome of the simulations 
at t = 0.25,0.75, 1.5 and 2.25tkh- We can directly see the 
strong impact these terms have on the stability of the sim- 
ulation. The standard a = 1.0, (3 = 2.0 is the most stable 
one, most probably due to the unphysical use of the a bulk 
viscosity. The use of a = and f3 = 2.0 or the Balsara 
switch renders very similar visual results. This is because 
the Balsara switch turns of viscosity where Vxt> is signifi- 
cant, which is the case for shearing flows like on the surface 
of the cloud. Note that Vxv is a very noisy quantity when 
measured using only 32 neighbours. By further lowering the 
shock capturing (3 viscosity we make the cloud even more 
unstable but it is not clear how physical this solution is. 
The shock front gets more blurred and we see strong post 
shock ringing effects. The reason for the increased instabil- 
ity in the a — 0, (3 = 0.5, and a = 0, /3 = 0.1 case is most 
probably due to high speed particles traveling through the 
poorly captured shock region and transferring momentum 
inside the cloud, perturbing it in an unphysical way. 

We see from these simulations how lowering viscosity 
will make the cloud less stable, which is expected from linear 
analysis. However, we still can't obtain agreement with the 
grid based codes. This leads us to suspect that there are 
more fundamental reasons behind the discrepancies. 

6.2 Resolving instabilities 

In order to create an even simpler test problem to compare 
instabilities between codes, we carried out a classical Kelvin- 
Helmholtz test using Gasoline and Enzo. We looked at the 
shearing motion of two gases of different densities and with 
small perturbations imprinted at the boundary. This cap- 
tures some of the hydrodynamics at the surface of the cloud 
in the blob test. 

The setup is a periodic box with dimensions 
{L x .L y , L z } = {1,1,1/32}, divided into two regions: one 
cold, high density and one warm, low density. The den- 
sity and temperature ratio is % = pb/pt — T t /Tb = c^/c^. 



putting the whole system in pressure equilibrium. The two 
layers are given constant and opposing shearing velocities, 
with the top layer moving leftward at a Mach number 
Mt = vt/ct « 0.11 and the bottom layer moving rightward 
at a Mach number Mb = Vb/cb ~ 0.34 in the case of \ — 10. 
The shear velocity becomes v s hear =0.68 Cb and the subsoni c 
regime will assure growth of instabilities (|Vietri et al.lll997F ) . 
This setup should mimic the growth of instabilities on the 
cloud surface. 

To trigger instabilities we have imposed sinusoidal per- 
turbation on the vertical velocity of the form 

v y (x) — Sv y sin(A27rx), (15) 

where Sv y is the amplitude of the perturbation in terms 
of the sound speed Cb and A is the wavelength of the mode 
which we have put to 1/6 in all of out tests. The perturbation 
is limited to a central strip around the interface of thickness 
5% of the box size. 

The initial conditions are again first generated using 
particles which are then mapped to the grid so that both 
codes have very similar starting points. An important issue 
for this type of test is how particles are distributed since this 
will introduce a certain amount of noise via discreteness. The 
most common techniques for this are: 

• Lattice: Particles ordered in a perfect grid. This min- 
imizes local density fluctuations. This type of IC is optimal 
for grid codes. 

• Poisson: Particles are randomly distributed which gen- 
erates local density variations, causing spurious pressure 
forces. 

• Glass: Particles are heated and relaxed until equilib- 
rium and homogeneity is found. This type of simulated an- 
nealing will create a relaxed system and is more optimal for 
SPH than for grids. 

Any initial condition that has local density variations will 
trigger small scale KH instabilities. We carried out this test 
using all three methods in order to illustrate their impact. 
The lattice is obviously perfect for grid codes, making a 
perfectly homogeneous gas. This quality does not automat- 
ically produce clean SPH initial conditions due to the aver- 
aging over nearby particles. The poisson ICs are very noisy 
in both the grid and SPH case, even though grid codes tend 
to smooth the noise over the cell sizes. The glass IC is in- 
tuitively the closest IC for both methods producing a self- 
consistent and homogeneous initial state for SPH simula- 
tions while leaving only small fluctuations for both grid and 
SPH methods. 

This set of simulations and their characteristics are 
summarized in Table Eland Fig. ll II shows the results, from 
top to bottom, GRID1, GRID3 and SPH3. We choose to 
show only one of the SPH results since all of these runs give 
the same result. GRID1 and GRID3 illustrate the difference 
between a highly idealised smooth setup (GRID1) and one 
with small scale noise (GRID3). 

GRID1 nicely produces the KH instabilities and the 
growth time is in excellent agreement with that expected 
from Eq.|H] This growth is not as clean in GRID3, which is 
to be expected due to local noise in density which alters the 
visual outcome. However the KHI is still well resolved and 
the growth time is comparable to the analytical expectation. 

The outcome of the SPH simulation is again very dif- 
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ferent from the grids. Perturbations are damped out very 
quickly both in velocity and density regardless of the initial 
conditions, the resolution and the viscosity. We conclude 
that SPH in the form used in astrophysical simulations to- 
date is unable to capture dynamical instabilities such as KH 
when density gradients are present. As we will show in the 
next section, the reason for this stems from the way hy- 
drodynamical forces are calculated in SPH in regions with 
strong gradients. 

6.3 Mind the gap 

Fig. 1121 shows a closeup of the SPH particles at the in- 
terface of the two fluids in SPH3 at t = tkh- There 
is a gap between them that has the size of the SPH 
smoothing kernel. This gap repeats periodically in each 
fluid, being smaller in the higher density fluid since the 
smoothing length (mean distance to the nearest 32 parti- 
cles) is smaller there. This feature is found in all of our 
SPH KH simulations. It occurs very quickly and becomes 
more prominent with time. This phen omenon has been 
discussed before in the literature (e.g. iRitchie &: Thomas! 
120011 : iTittlev et all feOOlft. especially in the context of nu- 
merical overcooling (IPearc e et al.lfl999HTh acker et all boOO ; 
ISoringel fc Hernauistl2002l : lMarri fc WhitJ200fl but no re- 
lation to resolving instabilities has been mentioned. 

The gap can also be clearly seen in the cloud test sim- 
ulation Fig. |5] Even though the gas is streaming with high 
velocity onto the leading surface of the cloud, the spurious 
pressure forces prevent it from making any physical contact. 
The reason that the cloud loses mass in the SPH simulation 
is due to the vacuum behind the cloud into which the cloud 
expands from its edges. Here the gradients become smooth 
and the gas can be removed by the pressure difference be- 
tween the cloud and the ambient medium that streams past. 

The effect can be explained in the following way: Eq. 1101 
is the force on each SPH particle coming from the summa- 
tion over the 32 nearest neighbours. The pressure is given 
by P rsj pT in the assumed case of an ideal gas. This force 
calculation formally assumes that temperature, and more 
importantly, density gradients are small within the smooth- 
ing kernel, where temperature is a quantity accumulated 
over time while density usually is re-estimated at each time 
step. When a particle from a hot low density region ap- 
proaches a cold high density region it will suddenly find a 
lot of neighbours at the edge of the smoothing sphere within 
the dense medium and its density will be overestimated. 
This leads to, through momentum conservation, a repulsive, 
fictitious, force on the particle, causing it to bounce back 
into the low density region. This behaviour leads to the for- 
mation of a gap between the two phases of the size ryj *2ihij , 
where hij is the effective smoothing kernel length, either ob- 
tained by using smoothin g length or smoothing kernel aver- 
aging (jHernauist &; Katzlll989r) , depending on the SPH im- 
plementation. Hot particles close to this gap will now have 
a strongly assymetric distribution of particles around them 
resulting in an average pressure force pointing back into the 
vacuum layer. Particles then travel back into the empty re- 
gion and the whole process is repeated. 

As mentioned above, this erroneous treatment of den- 
sity contrasts has also been fou nd to produce over cooling: in 
galaxy formation simulations. ITittlev et alJ (|200lh showed 



that in subsonic regimes this behaviour leads to fictitious ac- 
cretion of particles on the lateral sides of gas clouds such as 
the simulations showed in this paper. Solutions to this prob- 
lem has been attempted by several authors (see references 
above) by reformulating SPH to more accurately treat the 
particle interactions at steep boundaries. While this seems 
to remove the vacuum layer to some extent, it is unclear how 
this will affect the simulations discussed here. Possible so- 
lutions to the gap problem, such as modifying the viscosity 
weighting kernel, wil l be presented in a follow up paper by 
lAgertz et. all fcOOfl . 

That the gap is the origin of instability supression be- 
comes even more aparent by studying the KHI using a den- 
sity contrast x — 1? i n which the gap can not form. With this 
vanishing density gradient, SPH is able to capture the KHI, 
see Fig. [13] The left panel shows the KHI at t = tkh for the 
standard a = 1.0, f3 = 2.0 setting and the right panel shows 
the same timestep but using a — 0.01 and f3 = 1.0. The less 
evolved standard viscosity simulation points out the effects 
of viscosity discussed in section ^TTI Similar results have been 
recently found bv lJunk et. al. I (|20Q6"). 



7 SUMMARY 

In this paper we have carried out hydrodynamical simula- 
tions of a cold gas cloud interacting with an ambient hot 
moving gas using state of the art simulations codes. Strik- 
ing differences were found between the two main techniques 
for simulating fluids. While grid codes are able to resolve and 
treat dynamical instabilities and mixing, these processes are 
poorly or not at all resolved by the current SPH techniques. 
We show that the reason for this is that SPH, at least in 
the standard usage and formulation, inaccurately handles 
situations where density gradients are present. In these sit- 
uations, SPH particles of low density close to high density 
regions suffer erroneous pressure forces due to the assymet- 
ric density within the smoothing kernel. This causes a gap 
between regions of high density contrast. 

This behaviour has implications for many astrophysical 
situations. For example the stripping of gas from galaxies 
moving through a gaseous medium has already been dis- 
cussed in the literature. The origin of disk galaxies is an im- 
portant unsolved problem. Perhaps the inability to disrupt 
accreting gas clouds is one reason why numerical calcula- 
tions have failed to produce pure disk systems. Simulating 
star formation regions and feedback processes also relies on 
the correct ability to model turbulence and interacting mul- 
tiphase fluids. 

It should be noted that the behaviour of the grid and 
SPH mehods agree on timescales shorter than those of typi- 
cal dynamical instabilities such as the Kelvin-Helmholtz and 
Ray leigh- Taylor instabilities. In our specific test of a cold 
cloud engulfed in a hot wind, there is good agreement in 
the early gas stripping phase occuring due to pressure dif- 
ferences arising in the Bernoulli zones. As soon as the large 
scale instabilities have grown, the results of the different 
methods diverge. There are several possible solutions to this 
behaviour in SPH calculations which we will explore in a 
seperate work. 
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Figure 11. Density slices of, from top to bottom, GRID1, GRID3 and SPH3. The panels show the KH simulation at t = tkh/3, 2tkh/3 
and tkh • It is clear from the tests that the SPH method can not resolve the KH instability. 




Figure 12. A close up view of the SPH particles at the boundaries between the shearing layers (left) and closer zoom in (right) for 
SPH3 at tkh- We can clearly see artificial vacuum layers formed through erroneous pressure forces due to improper density calculations 
at density gradients. Even though the two fluids are moving relative to each other, the gap is so large that the viscosity information is 
not transferred between the two fluids. 



ACKNOWLEDGEMENTS 

We acknowledge support from the European Science Foun- 
dation who funded an exploratory workshop in Wengen 2004 
at which these tests were first discussed. FM and LM ac- 
knowledge support by the Swiss Institute of Technology 
through a Zwicky Prize Fellowship. OA would like to thank 
Alessandro Romeo and Peter Englmaier for valuable discus- 
sions. AJG acknowledge the support from the Polish Min- 
istry of Science through the grant 1P03D02626 and from 
the European Community's Human Potential Programme 
through the contract HPRN-CT-2002-00308, PLANETS. 
The AMR software (FLASH) used in this work was in part 
developed by the DOE-supported ASC / Alliance Centre for 
Astrophysical Thermonuclear Flashes at the University of 
Chicago. The FLASH calculations were performed at the In- 



terdisciplinary Centre for Mathematical and Computational 
Modeling in Warsaw, Poland. 



REFERENCES 

Abadi M. C, Moore B., Bower R. C, 1999, MNRAS, 308, 
947 

Agertz et. al. 2007, In preparation 

Balsara D. S., 1995, Journal of Computational Physics, 121, 
357 

Bryan G. L., Norman M. L., 1997, in Clarke D. A., West 
M. J., eds, ASP Conf. Ser. 123: Computational Astro- 
physics; 12th Kingston Meeting on Theoretical Astro- 
physics Simulating X-Ray Clusters with Adaptive Mesh 
Refinement, pp 363 — h 



Simulating fluids using SPH and grid techniques 15 




Figure 13. A zoom in of the SPH particles at the boundaries between the shearing layers for the isodensity SPH run with standard 
viscosity (left) and low viscosity (right) at tkh • The black and white region are particles that belonged to the initially separated shearing 
layers. We clearly see that the growth of the KHI in the standard implementation of SPH, and even stronger for the low viscosity version. 



Chandrasekhar S., 1961, Hydrodynamic and hydromag- 
netic stability. International Series of Monographs on 
Physics, Oxford: Clarendon, 1961 
Colella P., Glaz H. M., 1985, JCPh, 59, 264 
Colella P., Woodward P. R., 1984, JCPh, 54, 174 
Dolag K., Vazza F., Brunetti C, Tormen C, 2005, MN- 
RAS, 364, 753 

Doroshkevich A. C, Zeldovich I. B., 1981, Zhurnal Eksper- 
imental noi i Teoreticheskoi Fiziki, 80, 801 

Faber T. E., 1995, Fluid dynamics for physicists. Cam- 
bridge, New York: Cambridge University Press, — cl995 

Fryxell B., Olson K., Ricker P., Timmes F. X., Zingale M., 
Lamb D. Q., MacNeice P., Rosner R., Truran J. W., Tufo 
H., 2000, ApJS, 131, 273 

Gregori C, Miniati F., Ryu D., Jones T. W., 2000, ApJ, 
543, 775 

Hernquist L., Katz N., 1989, ApJS, 70, 419 
Inogamov N. A., 1999, Astrophysics and Space Physics Re- 
views, 10, 1 

Jones T. W., Ryu D., Tregillis I. L., 1996, ApJ, 473, 365 
Junk et. al. 2006, In preparation 
Khokhlov A. M., 1998, JCPh, 143, 519 
Klein R. I., McKee C. F., Colella P., 1994, ApJ, 420, 213 
Kravtsov A. V., 1999, PhD thesis, New Mexico State Uni- 
versity 

Kravtsov A. V., Klypin A., Hoffman Y., 2002, ApJ, 571, 
563 

Laney C, 1998, Computational gasdynamics. CUP, Cam- 
bridge 

Lattanzio J. C, Monaghan J. J., Pongracic H., Schwarz 
M. P., 1986, j-SIAM-J-SCI-STAT-COMP, 7, 591 

Leveque R. J., 1998, in Steiner O., Gautschy A., eds, Saas- 
Fee Advanced Course 27: Computational Methods for As- 
trophysical Fluid Flow. Nonlinear Conservation Laws and 
Finite Volume Methods, p. 1 

Marri S., White S. D. M., 2003, MNRAS, 345, 561 

Miniati F., Colella P. J., 2006, Comp. Phys. submitted, 
e-print:astro-ph/0608156 

Miniati F., Jones T. W., Ryu D., 1999, ApJ, 517, 242 

Monaghan J. J., 1992, ARA& A, 30, 543 

Mori M., Burkert A., 2000, ApJ, 538, 559 

Murray S. D., White S. D. M., Blondin J. M., Lin D. N. C, 
1993, ApJ, 407, 588 

Nulsen P. E. J., 1982, MNRAS, 198, 1007 

Pearce F. R., Jenkins A., Frenk C. S., Colberg J. M., White 
S. D. M., Thomas P. A., Couchman H. M. P., Peacock 
J. A., Efstathiou C, The Virgo Consortium 1999, ApJL, 
521, L99 

Quilis V., Moore B., Bower R., 2000, Science, 288, 1617 
Ritchie B. W., Thomas P. A., 2001, MNRAS, 323, 743 
Shu F. H., 1992, Physics of Astrophysics, Vol. II. Physics of 



Astrophysics, Vol. II, by Frank H. Shu. Published by Uni- 
versity Science Books, ISBN 0-935702-65-2, 476pp, 1992. 
Springel V., 2005, MNRAS, 364, 1105 
Springel V., Hernquist L., 2002, MNRAS, 333, 649 
Springel V., Yoshida N., White S. D. M., 2001, New As- 
tronomy, 6, 79 
Stadel J. C, 2001, Ph.D. Thesis 

Thacker R. J., Tittley E. R., Pearce F. R., Couchman 
H. M. P., Thomas P. A., 2000, MNRAS, 319, 619 

Tittley E. R., Pearce F. R., Couchman H. M. P., 2001, ApJ, 
561, 69 

van Leer B., 1979, JCPh, 32, 101 

Vietri M., Ferrara A., Miniati F., 1997, ApJ, 483, 262 
Vikhlinin A., Markevitch M., Murray S. S., 2001, ApJL, 
549, L47 

Wadsley J. W., Stadel J., Quinn T., 2004, New Astronomy, 
9, 137 

Watkins S. J., Bhattal A. S., Francis N., Turner J. A., 
Whitworth A. P., 1996, AApS, 119, 177 



